Does the Chapman— Enskog expansion for sheared granular gases converge? 
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The fundamental question addressed in this paper is whether the partial Chapman-Enskog ex- 



pansion Pxy 



X^fcLo '7'= {dux/dy) of the shear stress converges or not for a gas of inelastic 



hard spheres. By using a simple kinetic model it is shown that, in contrast to the elastic case, the 
above series does converge, the radius of convergence increasing with inelasticity. It is argued that 
this paradoxical conclusion is not an artifact of the kinetic model and can be understood in terms 
of the time evolution of the scaled shear rate in the uniform shear flow. 



PACS numbers: 45.70.Mg, 47.50.-d, 05.20.Dd 51.10.+y 

The hydrodynamic description of conventional fluids is 
usually restricted to the Navier-Stokes (NS) constitutive 
equations [i|. For instance, if the flow is incompressible 
(i.e., V • u = 0, where u is the flow velocity), Newton's 
law establishes a linear relationship P^y = —ijodux/dy 
between the shear stress Pxy and the shear rate dux/dy, 
where 770 is the shear viscosity and it has been assumed 
that duy/dx = 0. The NS constitutive equations rep- 
resent excellent approximations in most of the physical 
situations of experimental interest, even if the regime is 
turbulent [4|. In the case of a dilute gas, they can be 
justified under the assumption that the smallest of the 
characteristic hydrodynamic lengths (L) associated with 
the gradients of density (n), temperature (T), and fiow 
velocity (u) is much larger than the mean free path I 
of the gas particles, i.e., /j, = i/L <C 1. As a matter 
of fact, the Chapman-Enskog (CE) method provides a 
systematic scheme to obtain the normal solution of the 
Boltzmann equation as an expansion in powers of the 
uniformity parameter (or Knudsen number) fi [3[. The 
leading terms in the CE expansion yield the NS consti- 
tutive equations, with the bonus of providing expressions 
for the transport coefficients (like the shear viscosity 770) 
in terms of the microscopic properties and of the hydro- 
dynamic quantities. 

A fundamental question concerning the CE method is 
the nature (convergent versus divergent) of the expan- 
sion in powers of ^. More than 40 years ago Grad Q 
provided compelling arguments on the asymptotic char- 
acter of the CE expansion. This means that, as expected 
on physical grounds, the CE series truncated at a given 
order (e.g., NS, Burnett, super-Burnett, ...) becomes 
closer and closer to the true value as /i becomes smaller 
and smaller. Of course, the asymptotic character of the 
CE expansion does not imply (but is implied by) the 
stronger condition of convergence. However, McLennan 
[5| was able to prove convergence of the partial sum of the 
CE series made of linear terms (i.e., terms of the form 



V,,V,, 



, Vi^A, where A = n, T, or u) for a general 
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class of cutoff potentials. 

What about the nonlinear terms of the series? To be 
more specific, let us consider the following subclass of the 
full CE series of the shear stress: 



Pxv = -^Vk {dux/dyf^^^ 



xy 



fc=0 



(1) 



where ryo is the NS shear viscosity, rji is a super-Burnett 
coefficient, and so on. The full CE series of Pxy reduces 
to the partial series ([T]) if (and only if) dui/dxj = jSixSjy 
and Vn = VT = V7 = 0, i.e., the only nonzero hydro- 
dynamic gradient is a uniform shear rate 7 = dux/dy. 
Interestingly enough, there exists a physical state (the 
so-called simple or uniform shear flow, USE) that is 
consistent with those conditions [a, Q- In that state, 
the density and the shear rate are constant in time 
but the temperature increases due to viscous heating. 
The identification of the characteristic hydrodynamic 
length L is (except for a numerical factor) unambiguous: 
L ^ •y/2T/?Ti/7 (where m is the mass of a particle), so 
that the uniformity parameter becomes /i = i/v, where 
v ^ ^J2T /m/ 1 is a characteristic collision frequency. 
Thus, the CE expansion ([T|) can be recast into the di- 
mensionless form 



Pxy/nT=-fiF{p^), F(z) 



k=Q 



(2) 



where Ck = {rjk / nT)v'^''^^ . Despite its simple definition, 
the USE state is complex enough to prevent an exact so- 
lution of the nonlinear Boltzmann equation. However, 
the problem becomes solvable in the framework of the 
Bhatnagar-Gross-Krook (BGK) model kinetic equation 
[3] and the solution shows that, for a wide class of re- 
pulsive potentials, the CE expansion ([T]) is divergent [^. 
More specifically, the dimensionless coefficients c^ behave 
for large k as \ck\ ~ {2/dYk\ in the case of d-dimensional 
hard spheres. 

In the preceding paragraphs we have considered con- 
ventional gases made of particles that collide elastically. 
On the other hand, the same issues discussed above can 
be applied to granular gases, i.e., large assemblies of 
(mesoscopic or macroscopic) particles which collide in- 
elastically and are maintained in fluidized states. Apart 



from their practical interest, granular gases are physi- 
cal systems worth studying at a fundamental level be- 
cause they are intrinsically out of equilibrium and ex- 
hibit a wide repertoire of complex and exotic behavior 
[9|, |lO|. Nonequilibrium statistical-mechanical concepts 
and tools, in particular the kinetic theory approach based 
on the Boltzmann and Enskog equations suitably modi- 
fied to account for inelastic collisions, have proven to be 
useful to understand the behavior of granular gases [ll| . 
Most of the theoretical efforts have focused on a minimal 
model of granular gases consisting of smooth inelastic 
hard spheres characterized by a constant coefficient of 
normal restitution a < 1. Specifically, the CE method 
has been applied to the (inelastic) Boltzmann and En- 
skog kinetic equations and the NS transport coefficients 
have been derived 12]. 

The question I want to address in this Letter is, does 
the nonlinear subclass of the full CE expansion converge 
in the case of a gas of inelastic hard spheres? Given that 
the answer is negative when the gas is made of elastic 
hard spheres (a = 1) [a], the strong challenges to the 
validity of hydrodynamics in granular media [lO|, and 
the inherently non-Newtonian nature of the steady USE 
of granular gases [IjI, it seems plausible to expect that 
the (partial) CE series ^ is divergent for granular gases. 
It will be shown below that, paradoxically, the series ^ 
does converge in the case of inelastic hard spheres and 
that, in fact, the radius of convergence increases with 
increasing inelasticity. 

The energy balance equation for inelastic hard spheres 
in the USE is 



five collision frequency for elastic spheres. The cooling 
term on the right-hand side of Eq. ([3]) competes with the 
viscous heating term, so that, depending on the initial 
state and the value of 7, the temperature either grows 
or decreases with time until a steady state is eventually 
reached [ij, llJ] . In order to relate the shear stress Pxy (t) 
to the shear rate 7 and to the temperature T{t), and thus 
analyze the series representation ©, let us replace the 
Boltzmann equation by the BGK-like kinetic model [1^ 

{dt + v V)/ = ^u{f - /o) + (C/2)av • (V/), (4) 

where / is the velocity distribution function, V = v — u 
is the peculiar velocity, /o is the local version of the ho- 
mogeneous cooling state distribution [ll| (parameterized 
by the actual fields n, u, and T), and v is an effective 
collision frequency. Here the simple choice v 



l+a ,, ,„ 
-^-1^0 IS 

made, so that e = (^/v = ^^(1 ~ ct)- Taking moments 
in Eq. (|4]) one gets 

dtP.yit)^ -^Pyyit)-Mt)+at)]P,yit), 

dtPyyit) = ni^{t)T{t) - [u{t) + at)]Pyy{t),^ ^ ' 

where Pyy is a normal stress. Note that the explicit form 
of /o is not needed in the derivation of Eq. ([5]) and so 
no assumption of / or /o being close to a Maxwellian is 
taken. Equation (J5]), with v = i±^ [l - ^(1 - a)] uq, 
is also obtained from the Boltzmann equation in Grad's 
moment approximation [l3|. Equations ([3]) and (O con- 
stitute a closed set of equations for the evolution of T, 
Pxy^ and Pyy In order to describe the hydrodynamic 
regime and analyze the CE expansion we must focus on 
the nonlinear dependence of the scaled viscosity function 



dtT{t) = -{2/dn)jPxy{t) - ({t)T{t), (3) F{^^^it)) = -Pxy{t)/nT{t)n{t) as a function of the scaled 



where (^ is the so-called cooling rate jll| . Erom the Boltz- 
mann equation it is possible to show that it is approx- 



imately given by C 



d+2 
id 



(1 



)va; Vf) being an cffec- 



shear rate (or uniformity parameter) /i(t) = ^jvit). 
Elinfination of time in favor of z = /i^ oc l/T in Eqs. 
([3]) and (O yields the following single second-order ordi- 
nary differential equation: 



l-(l + ^.F)^^4(e 



lzF){zdA{^-\zF)[F- 



2zd,F)\ + (\ 



V 



-^zF){F- 



2zd^F) + 2zd^[{l 



^zF)F]}. (6) 



Insertion of the series expansion F{z) ~ X]fc°=o ^kz'^ into 
Eq. ([5]) allows one to get a recursive relation expressing 
the coefficient Ck for fc > 1 in terms of the lower-order 
coefficients {ck',0 < k' < fc — 1}. With the help of a 
computer algebra system, one can obtain the coefficients 
Ck, for given d and a, as exact rational numbers up to 
a maximum value of fc limited by computer time and 
internal memory. The results of Ck presented here have 
been obtained for k < 400. 

The first few coefficients Ck are listed in Table [l] for 
d = 3 and a = 0.5, 0.7, 0.9, 0.99, and 1 (elastic hard 
spheres). As is already apparent from Table U the mag- 
nitude of the coefficients increases as the system becomes 



less inelastic, this effect being more and more dramatic 
with increasing fc. Eor instance, the value of |c2o| is 
1.2 X 10-^ 5.2 X 10-3, 2.1 X 10^ 1.4 x 10", and 1.1 x lO^^ 
for a = 0.5, 0.7, 0.9, 0.99, and 1, respectively. Figure 
[T] shows a log- normal plot of \ck\ for the four inelastic 
systems considered. The results are consistent with an 
asymptotic behavior of the form \ck\ ~ Ak~°'z~'', where 
Zs = |e(l -I- e)^ is the steady-state value of the square 
uniformity parameter. The exponent a and the ampli- 
tude A are estimated by a linear fit of a plot of ln(|cfc|Zj ) 
versus In fc for 100 < fc < 400, the fit being less robust in 
the case a = 0.5 than in the other cases. The numerical 
values of Zs , as well as of the estimates of a and In A are 



TABLE I: CE coefficients c^ for < fc < 6, steady-state value of the square uniformity parameter, Zs, and estimates of the 
parameters a and In^l associated with the asymptotic behavior of |cfc|. 



Co 



Cl 



C2 



C3 



C4 



C5 



C6 



In A 



0.5 
0.7 
0.9 
0.99 
1 



0.827586 
0.888889 
0.960000 
0.995851 
1 



-0.380155 
-0.536326 
-0.798195 
-0.976777 

-1 



0.0677085 
0.146688 
0.377471 
0.627438 
0.666667 



0.0404196 
0.110743 
0.376385 
0.720350 

0.777778 



-0.00251705 

-0.0180970 

-0.188551 

-0.694311 

-0.814815 



-0.0119767 

-0.0727681 

-0.778076 

-3.20037 

-3.82716 



-0.00399593 

-0.0336082 

-0.618849 

-3.75620 

-4.72428 



1.25434 

0.585938 

0.146701 

0.0127092 
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FIG. 1: (color online) Plot of sgn(cfc) In |cfc| for d = 3 and 
a — 0.5, 0.7, 0.9, and 0.99. The thin solid lines represent 
±(— fclnzs — alnfc -I- In j4), where a and InA are fitting pa- 
rameters. 




FIG. 2: (color online) Scaled viscosity function F{z) for d — 3 
and a = 0.5. The thin solid line represents the numerical 
solution of Eq. (|6]), while the dashed lines represent the trun- 
cated CE expansion Fn{z) with N = 5, 20, 50, 100, and 
200. The thick solid line represents the results obtained from 
Monte Carlo simulations of the Boltzmann equation [I^. The 
steady-state points {zs,Fs) obtained from the kinetic model 
and the simulations are indicated by the triangle and the cir- 
cle, respectively. The inset includes Fn{z) with A'^ = 1-4 in 
the region < z < 1.5. 



also included in Table HI In principle, the viscosity func- 
tion F{z) is different for each dimensionality. However, 
inspection of Eq. ([6]) shows that, for a given value of e, 
the dependence of F on z and d occurs only through the 
combination z/d. Therefore, seen as a function of e and 
z/zs, F{z) is independent of d. 



The large-fc asymptotic behavior ln|cfe| '^ — /clnzs for 
a < 1 implies that the CE series ([T]) is convergent for 
inelastic hard spheres. Moreover, the radius of conver- 
gence is dux/dy = vzs and thus increases with inelas- 
ticity. Equation ^ shows that z = z^ is indeed a singular 
point since the steady state {zs,Fs), with Fg = (1 -|-e)^^, 
is the solution of 1 - (l + |zF)^ F = e-^zF = 0. What 
is relevant here is that z = is a regular point (except 
in the elastic case) and z = Zs is the singularity of F(z) 
in the complex z plane closest to the origin. The con- 
vergent character of the series ([2]) for the most inelastic 
case considered here (a = 0.5) is illustrated in Fig. [21 
where the numerical solution [13l| of Eq. © is compared 

with the truncated CE series Fn{z) = X]fc=o '^k^'' for sev- 
eral values of the truncation order N. One can observe 
that the truncated series agree among themselves and 
with the numerical solution for z < Zg. For z > z^, how- 
ever, the CE expansion becomes useless and one must de- 
termine F{z) numerically from Eq. ^. An alternative 
method, successfully used in the elastic case Q, would 
consist of expanding F{z) around the point at infinity as 
F{z) = z-^^^J2kLoCkZ-''/^. Figure El also includes re- 
sults recently obtained from Monte Carlo simulations of 
the Boltzmann equation ^14'|, which show that the pre- 
dictions of Eq. ^ are quantitatively accurate, even for 
strong inelasticity, except that Eq. ^ underestimates the 
location of the steady-state point {zs,Fs). The simula- 
tion curve shown in Fig.[2]is the collapse of data obtained 
by the direct simulation Monte Carlo method starting 
from 20 different initial conditions (ten for z < Zs and 
ten for z > Zs), letting the system evolve in time, and 
discarding the kinetic transients lasting a few mean free 
times. This explains why the simulation curve does not 
reach the point z = 0. For simulation details the reader 
is referred to Ref. [ij. 

Is the paradoxical regularization by inelasticity of the 
CE series ([T|) an artifact of the kinetic model? The follow- 
ing physical argument suggests that this is not the case. 
Let us assume that the reference homogeneous state is 
perturbed at time i = by a weak shear rate 7. Thus the 
homogeneous state (/x — 0) is stable or unstable against 
this USE perturbation depending on whether the time- 
dependent uniformity parameter /x(t) = ^ /v{t) asymp- 
totically goes to zero or grows with time. The first situ- 
ation occurs in a conventional gas of elastic hard spheres 
since the viscous heating produces a monotonic increase 
of temperature [cf. Eq. ([3]) with C = 0] and thus /i(t) -^ 0. 



However, in a gas of inelastic hard spheres it is always 
possible that the perturbation is small enough to make 
the cooling rate prevail over the viscous heating rate and 
so the temperature keeps decreasing in time [cf. Eq. ([3]) 
with C > 0] until the steady state is eventually reached; 

1/2 

thus the uniformity parameter grows, i.e., ii{t) — > Zg . 
Since the CE expansion ^ measures the departure from 
the reference homogeneous state (/u = 0), it is reason- 
able to expect that the series diverges if /i(t) goes to zero 
while it converges if fi(t) increases with time. 

The above heuristic argument can also be applied to 
the simplest CE series in a compressible flow, namely 



k=0 



(7) 



The full CE expansion of Pyy reduces to the series 
([7]) in the uniform longitudinal flow (ULF) charac- 
terized by jtU?] duy/dy = 7(t) — 7o/(l + jot) and 
n(t) = {no/jo)j{t). The exact energy balance equa- 
tion, dtT{t) = -{2/dn)j{t)Pyy{t) ~ C{t)T{t), shows that 
the uniformity parameter /i(i) = j{t)/i'{t) increases with 
time if 70 > 0, both for elastic and inelastic collisions. 
However, if 70 < 0, |/x(i)| — > in the elastic case whereas 
/i(i)| increases toward a stationary value in the inelastic 
case. Thus, the CE expansion ([7]) is expected to diverge 
for conventional gases and converge for granular gases. 
This expectation is confirmed by an analysis of the ULF 
based again on the kinetic model ([4]). 

It is important to note that the dependence on n, T, 



and a of the transport coefficients rjk and ry^. in the partial 
CE expansions ((T|) and , respectively, is not influenced 
by the specific state under consideration. Accordingly, 
the series ([T]) and ([7]) converge or diverge irrespective of 
whether the system is in the USE, the ULF, or in any 
other state. The advantage of the USE and ULE is that 
the full CE series of the shear stress and the normal stress 
reduce to the partial series H]) and ([7]) , respectively, thus 
allowing us to explore their character in a rather detailed 
way. Regarding the CE subseries made of linear terms, it 
is reasonable to expect that, as in the elastic case, it also 
converges for inelastic collisions, especially if one takes 
into account the exact mapping between the inelastic and 
elastic versions of the Enskog-Lorentz equations [131 ■ 

The CE expansion is the main route to hydrodynamics 
and so its convergence or divergence has a renewed inter- 
est in granular gases in view of some debate on the appli- 
cability of^hydrodynamics to this class of nonequilibrium 
systems [10]. I expect that this Letter can contribute to 
a clarification of this controversial issue by presenting a 
case study where the application of the CE expansion to 
describe the nonlinear regime might have a larger prac- 
tical interest in granular than in conventional gases. 
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